1/f Noise Project Manuscript 1 Final Analysis
Fine-Scale Spatial Heterogeneity in Phytochemistry Increases Herbivore Performance by Enhancing Herbivore Foraging
Vincent S. Pan, Enakshi Ghosh, Paul J. Ode, William C. Wetzel, Kadeem J. Gilbert, Ian S. Pearse
May 29, 2024
1 Setup
1.1 Load libraries
library(survival)
library(piecewiseSEM)
source("spat1f/init_analysis.R")Data last updated on 2024-05-031.2 Handy functions not in spat1f
trend_by_pre_wt <- function(model, term){
emmeans::emtrends(model, var = "cat_pre_wt_log_scale", specs = term) %>%
pairs()
}
contrast_by_pre_wt <- function(model, term, type = c("pairwise", "trt.vs.ctrl"), at = c(-2, 2)){
type <- match.arg(type)
f <- as.formula(sprintf("%s ~ %s | cat_pre_wt_log_scale", type, term))
emmeans::emmeans(model,
f,
var = term,
at = list(cat_pre_wt_log_scale = at))
}
check_mod <- function(model){
DHARMa::simulateResiduals(model) %>%
plot()
performance::posterior_predictive_check(model)
}1.3 Quick cleaning
# Main data.frame for analyses
d <- ref_data %>%
filter(error == 0) %>%
mutate(
cat_pre_wt_log = log(cat_pre_wt),
cat_pre_wt_log_scale = as.numeric(scale(log(cat_pre_wt))),
cat_size = ifelse(cat_pre_wt < median(cat_pre_wt), "small", "big"),
has_var = ifelse(var_trt != "constant", 1, 0),
mean_toxic_conc_scale = as.numeric(scale(mean_toxic_conc)),
mean_toxic_scale = as.numeric(scale(mean_toxic)),
area_herb_log_scale = as.numeric(scale(log(area_herb+1))),
var_toxic_12_scale = as.numeric(scale(var_toxic_12)),
on_toxic_logit_scale = as.numeric(scale(adjust_prop(on_toxic,
trans = "emp",
nudge.method = "none",
na.action = "ignore"))),
ava_qual_scale = as.numeric(scale(ava_qual)),
ava_qual_logit_scale = as.numeric(scale(adjust_prop(ava_qual,
trans = "emp",
nudge.method = "none",
na.action = "ignore"))),
sl_mean_obs_log_scale = as.numeric(scale(log(sl_mean_obs))),
cat_pre_wt_log_scale_sq = as.numeric(scale(log(cat_pre_wt)^2)),
RGR_scale = as.numeric(scale(RGR)),
sl_mean_pred1 = shape_1 * scale_1,
sl_mean_pred2 = shape_2 * scale_2,
sl_mean_pred3 = shape_3 * scale_3,
sl_mean_pred4 = shape_4 * scale_4,
k1 = scale_1 / scale_2,
k2 = scale_3 / scale_4,
prop_explore_logit = qlogis(prop_explore),
prop_explore_logit_scale = as.numeric(scale(prop_explore_logit))
)
# Subset of data.frame for SEM
d2 <- d %>%
filter(
var_trt != "constant"
) %>%
filter(
!is.na(area_herb_log_scale) &
!is.na(var_toxic_12_scale) &
!is.na(mean_toxic_conc_scale) &
!is.na(sl_mean_obs) &
!is.na(prop_explore) &
!is.na(RGR)) %>%
mutate(
var_high = ifelse(var_trt == "high_var", 1, 0),
beta_red = ifelse(as.character(beta) == 5, 1, 0),
beta_white = ifelse(as.character(beta) == 0, 1, 0),
beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
)
# Subset of data.frame for more detailed analyses
d3 <- d %>%
filter(
var_trt != "constant"
) %>%
filter_at(vars(contains("shape_"), contains("scale_[0-9]"), contains("sl_mean_pred")),
all_vars(. > 0)) %>%
filter(sl_mean_pred2 < 2000) %>%
mutate(
var_high = ifelse(var_trt == "high_var", 1, 0),
beta_red = ifelse(as.character(beta) == 5, 1, 0),
beta_white = ifelse(as.character(beta) == 0, 1, 0),
beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
)2 Herbivore Performance
2.1 RGR
2.1.1 Model summary
rgr_m <- glmmTMB(
RGR ~
(var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # residual plot show significant non-linearity
(1|session_id),
data = d %>%
filter(
var_trt != "constant"
)
); summary(rgr_m) Family: gaussian ( identity )
Formula:
RGR ~ (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
(1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
-879.2 -849.4 450.6 -901.2 100
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.979e-06 0.001995
Residual 1.603e-05 0.004004
Number of obs: 111, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 1.6e-05
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.0162504 0.0013457 12.075 < 2e-16 ***
var_trtlow_var -0.0034872 0.0007740 -4.506 6.62e-06 ***
beta0 -0.0003065 0.0009664 -0.317 0.75112
beta5 -0.0013545 0.0009173 -1.477 0.13978
cat_pre_wt_log_scale -0.0048759 0.0011991 -4.066 4.78e-05 ***
I(cat_pre_wt_log_scale^2) -0.0021953 0.0006731 -3.262 0.00111 **
var_trtlow_var:cat_pre_wt_log_scale 0.0025873 0.0007897 3.276 0.00105 **
beta0:cat_pre_wt_log_scale 0.0011924 0.0009491 1.256 0.20897
beta5:cat_pre_wt_log_scale 0.0024731 0.0009710 2.547 0.01087 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(rgr_m, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: RGR
Chisq Df Pr(>Chisq)
(Intercept) 145.8141 1 < 2.2e-16 ***
var_trt 20.3008 1 6.617e-06 ***
beta 2.3963 2 0.301758
cat_pre_wt_log_scale 16.5348 1 4.777e-05 ***
I(cat_pre_wt_log_scale^2) 10.6374 1 0.001108 **
var_trt:cat_pre_wt_log_scale 10.7334 1 0.001052 **
beta:cat_pre_wt_log_scale 6.4936 2 0.038899 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.1.2 Post-hoc contrasts
trend_by_pre_wt(rgr_m, "var_trt") contrast estimate SE df t.ratio p.value
high_var - low_var -0.00259 0.00079 100 -3.276 0.0014
Results are averaged over the levels of: beta contrast_by_pre_wt(rgr_m, "var_trt")$emmeans
cat_pre_wt_log_scale = -2:
var_trt emmean SE df lower.CL upper.CL
high_var 0.014224 0.00355 100 0.007178 0.02127
low_var 0.005562 0.00312 100 -0.000635 0.01176
cat_pre_wt_log_scale = 2:
var_trt emmean SE df lower.CL upper.CL
high_var -0.000393 0.00208 100 -0.004512 0.00373
low_var 0.001295 0.00238 100 -0.003421 0.00601
Results are averaged over the levels of: beta
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
high_var - low_var 0.00866 0.00181 100 4.790 <.0001
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
high_var - low_var -0.00169 0.00171 100 -0.988 0.3256
Results are averaged over the levels of: beta trend_by_pre_wt(rgr_m, "beta") contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.00119 0.000949 100 -1.256 0.4232
(beta-5) - beta5 -0.00247 0.000971 100 -2.547 0.0329
beta0 - beta5 -0.00128 0.000940 100 -1.363 0.3644
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates contrast_by_pre_wt(rgr_m, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df lower.CL upper.CL
-5 0.012890 0.00382 100 0.005313 0.02047
0 0.010199 0.00333 100 0.003591 0.01681
5 0.006589 0.00319 100 0.000255 0.01292
cat_pre_wt_log_scale = 2:
beta emmean SE df lower.CL upper.CL
-5 -0.001439 0.00226 100 -0.005922 0.00304
0 0.000639 0.00231 100 -0.003940 0.00522
5 0.002153 0.00257 100 -0.002946 0.00725
Results are averaged over the levels of: var_trt
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 0.00269 0.00221 100 1.218 0.4456
(beta-5) - beta5 0.00630 0.00223 100 2.820 0.0158
beta0 - beta5 0.00361 0.00211 100 1.711 0.2062
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.00208 0.00205 100 -1.015 0.5689
(beta-5) - beta5 -0.00359 0.00206 100 -1.746 0.1935
beta0 - beta5 -0.00151 0.00211 100 -0.719 0.7529
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates 2.1.3 Misc stats
r2(rgr_m, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.555
Marginal R2: 0.4442.1.4 Check model
check_mod(rgr_m)
2.1.5 Plot
g1 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) +
geom_line(aes(color = beta), linewidth = 2) +
geom_point(
data = rgr_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = beta),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("fill","color")) +
labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
color = expression(beta), fill = expression(beta))
g2 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","var_trt")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = var_trt), alpha = 0.2) +
geom_line(aes(color = var_trt), linewidth = 2) +
geom_point(
data = rgr_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = var_trt),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
theme(legend.position = "top") +
scale_color_discrete(type = c("navy","skyblue"),
aesthetics = c("fill","color"),
label = c("High", "Low")) +
labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
color = "Variation", fill = "Variation")
ggarrange(g1, g2 + labs(y = ""))
2.2 Time to pupation
2.2.1 Model summary
pupt_m_full <- glmmTMB(
time_to_pupation ~
(var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = Gamma(link = "log"),
); summary(pupt_m_full) Family: Gamma ( log )
Formula: time_to_pupation ~ (var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
338.1 366.3 -158.0 316.1 85
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 2.839e-11 5.328e-06
Number of obs: 96, groups: session_id, 5
Dispersion estimate for Gamma family (sigma^2): 0.0254
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.180302 0.036853 59.16 < 2e-16 ***
var_trtlow_var 0.102145 0.033998 3.00 0.00266 **
beta0 -0.001845 0.040823 -0.05 0.96395
beta5 0.059543 0.041181 1.45 0.14821
cat_pre_wt_log_scale -0.428662 0.032643 -13.13 < 2e-16 ***
I(cat_pre_wt_log_scale^2) -0.081267 0.019069 -4.26 2.03e-05 ***
var_trtlow_var:cat_pre_wt_log_scale 0.026518 0.038876 0.68 0.49516
beta0:cat_pre_wt_log_scale -0.030363 0.043195 -0.70 0.48210
beta5:cat_pre_wt_log_scale -0.126801 0.047095 -2.69 0.00709 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(pupt_m_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: time_to_pupation
Chisq Df Pr(>Chisq)
(Intercept) 3500.2207 1 < 2.2e-16 ***
var_trt 9.0264 1 0.002661 **
beta 2.7351 2 0.254736
cat_pre_wt_log_scale 172.4448 1 < 2.2e-16 ***
I(cat_pre_wt_log_scale^2) 18.1614 1 2.029e-05 ***
var_trt:cat_pre_wt_log_scale 0.4653 1 0.495162
beta:cat_pre_wt_log_scale 7.5589 2 0.022836 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pupt_m <- glmmTMB(
time_to_pupation ~
var_trt + beta * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = Gamma(link = "log"),
); summary(pupt_m) Family: Gamma ( log )
Formula:
time_to_pupation ~ var_trt + beta * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
(1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
336.5 362.2 -158.3 316.5 86
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 2.159e-11 4.647e-06
Number of obs: 96, groups: session_id, 5
Dispersion estimate for Gamma family (sigma^2): 0.0255
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.180555 0.036937 59.03 < 2e-16 ***
var_trtlow_var 0.107447 0.033197 3.24 0.00121 **
beta0 -0.003578 0.040842 -0.09 0.93019
beta5 0.059536 0.041271 1.44 0.14915
cat_pre_wt_log_scale -0.420147 0.030221 -13.90 < 2e-16 ***
I(cat_pre_wt_log_scale^2) -0.083365 0.018881 -4.42 1.01e-05 ***
beta0:cat_pre_wt_log_scale -0.026011 0.042849 -0.61 0.54383
beta5:cat_pre_wt_log_scale -0.126834 0.047183 -2.69 0.00719 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(pupt_m, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: time_to_pupation
Chisq Df Pr(>Chisq)
(Intercept) 3485.0791 1 < 2.2e-16 ***
var_trt 10.4759 1 0.001209 **
beta 2.8053 2 0.245947
cat_pre_wt_log_scale 193.2743 1 < 2.2e-16 ***
I(cat_pre_wt_log_scale^2) 19.4946 1 1.009e-05 ***
beta:cat_pre_wt_log_scale 7.6890 2 0.021397 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.2.2 Post-hoc contrasts
contrast_by_pre_wt(pupt_m, "var_trt")$emmeans
cat_pre_wt_log_scale = -2:
var_trt emmean SE df asymp.LCL asymp.UCL
high_var 2.808 0.0769 Inf 2.657 2.96
low_var 2.915 0.0818 Inf 2.755 3.08
cat_pre_wt_log_scale = 2:
var_trt emmean SE df asymp.LCL asymp.UCL
high_var 0.924 0.0672 Inf 0.792 1.06
low_var 1.031 0.0725 Inf 0.889 1.17
Results are averaged over the levels of: beta
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df z.ratio p.value
high_var - low_var -0.107 0.0332 Inf -3.237 0.0012
cat_pre_wt_log_scale = 2:
contrast estimate SE df z.ratio p.value
high_var - low_var -0.107 0.0332 Inf -3.237 0.0012
Results are averaged over the levels of: beta
Results are given on the log (not the response) scale. trend_by_pre_wt(pupt_m, "beta") contrast estimate SE df z.ratio p.value
(beta-5) - beta0 0.026 0.0428 Inf 0.607 0.8163
(beta-5) - beta5 0.127 0.0472 Inf 2.688 0.0197
beta0 - beta5 0.101 0.0473 Inf 2.132 0.0834
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates contrast_by_pre_wt(pupt_m, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df asymp.LCL asymp.UCL
-5 2.741 0.0941 Inf 2.557 2.93
0 2.790 0.0971 Inf 2.599 2.98
5 3.054 0.1088 Inf 2.841 3.27
cat_pre_wt_log_scale = 2:
beta emmean SE df asymp.LCL asymp.UCL
-5 1.061 0.0845 Inf 0.895 1.23
0 1.005 0.0805 Inf 0.847 1.16
5 0.866 0.0917 Inf 0.687 1.05
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df z.ratio p.value
(beta-5) - beta0 -0.0484 0.1018 Inf -0.476 0.8828
(beta-5) - beta5 -0.3132 0.1133 Inf -2.763 0.0158
beta0 - beta5 -0.2648 0.1136 Inf -2.331 0.0516
cat_pre_wt_log_scale = 2:
contrast estimate SE df z.ratio p.value
(beta-5) - beta0 0.0556 0.0875 Inf 0.636 0.8005
(beta-5) - beta5 0.1941 0.0915 Inf 2.122 0.0854
beta0 - beta5 0.1385 0.0928 Inf 1.492 0.2945
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates 2.2.3 Misc stats
r2(pupt_m, tolerance = 10^-10)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.8812.2.4 Check model
check_mod(pupt_m)
2.2.5 Plot
g3 <- marginal_effects(pupt_m, terms = c("cat_pre_wt_log_scale","beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) +
geom_line(aes(color = beta), linewidth = 2) +
geom_point(
data = pupt_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = time_to_pupation, color = beta),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log2") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("fill","color")) +
labs(x = "Cat pre-weight (g)", y = "Time to pupation (days)",
color = expression(beta), fill = expression(beta))
g4 <- marginal_effects(pupt_m, terms = c("var_trt")) %>%
ggplot(aes(x = var_trt, y = yhat)) +
geom_point(
data = pupt_m$frame,
aes(x = var_trt, y = time_to_pupation, color = var_trt),
position = position_jitter(width = 0.2),
size = 3, alpha = 0.8
) +
geom_pointrange(
aes(ymax = upper, ymin = lower),
color = "black",
size = 1, linewidth = 2, shape = 1,
) +
theme_bw(base_size = 15) +
scale_y_continuous(trans = "log2") +
theme(legend.position = "top") +
scale_color_discrete(type = c("navy","skyblue"),
aesthetics = c("color"),
label = c("High", "Low")) +
labs(x = "Variation", y = "Time to pupation (days)",
color = "Variation")
ggarrange(g3, g4 + labs(y = ""))
2.3 Eclosure
2.3.1 Model summary
eclosure_m_full <- glmmTMB(
eclosed ~
(var_trt + beta) * cat_pre_wt_log_scale +
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = binomial(link = "logit"),
); summary(eclosure_m_full) Family: binomial ( logit )
Formula:
eclosed ~ (var_trt + beta) * cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
185.3 210.9 -83.6 167.3 119
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.1661 0.4076
Number of obs: 128, groups: session_id, 5
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.53530 0.41050 1.304 0.192
var_trtlow_var 0.01707 0.37507 0.046 0.964
beta0 -0.21544 0.46265 -0.466 0.641
beta5 -0.50749 0.45596 -1.113 0.266
cat_pre_wt_log_scale 0.10452 0.36089 0.290 0.772
var_trtlow_var:cat_pre_wt_log_scale 0.17526 0.36489 0.480 0.631
beta0:cat_pre_wt_log_scale 0.24778 0.43553 0.569 0.569
beta5:cat_pre_wt_log_scale 0.20159 0.44699 0.451 0.652car::Anova(eclosure_m_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: eclosed
Chisq Df Pr(>Chisq)
(Intercept) 1.7004 1 0.1922
var_trt 0.0021 1 0.9637
beta 1.2502 2 0.5352
cat_pre_wt_log_scale 0.0839 1 0.7721
var_trt:cat_pre_wt_log_scale 0.2307 1 0.6310
beta:cat_pre_wt_log_scale 0.3659 2 0.8328eclosure_m <- glmmTMB(
eclosed ~
(var_trt + beta) + cat_pre_wt_log_scale +
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = binomial(link = "logit"),
); summary(eclosure_m) Family: binomial ( logit )
Formula:
eclosed ~ (var_trt + beta) + cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
179.9 197.0 -84.0 167.9 122
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.1676 0.4094
Number of obs: 128, groups: session_id, 5
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.551286 0.413578 1.333 0.183
var_trtlow_var 0.005064 0.373402 0.014 0.989
beta0 -0.237692 0.460689 -0.516 0.606
beta5 -0.523521 0.456587 -1.147 0.252
cat_pre_wt_log_scale 0.330169 0.220049 1.500 0.134car::Anova(eclosure_m, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: eclosed
Chisq Df Pr(>Chisq)
var_trt 0.0002 1 0.9892
beta 1.3212 2 0.5165
cat_pre_wt_log_scale 2.2513 1 0.13352.3.2 Post-hoc contrasts
2.3.3 Misc stats
r2(eclosure_m, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.092
Marginal R2: 0.0462.3.4 Check model
check_mod(eclosure_m)
2.3.5 Plot
2.4 Survival
2.4.1 Model summary
surv_m_full <- coxph(
Surv(surv_time, observed_dead) ~
(var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
strata(session_id),
data = d %>%
filter(var_trt != "constant") %>%
mutate(
beta = as.factor(as.character(beta))
),
); summary(surv_m_full)Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) *
cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id),
data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))
n= 128, number of events= 105
coef exp(coef) se(coef) z Pr(>|z|)
var_trtlow_var -0.1115 0.8945 0.2124 -0.525 0.599727
beta0 -0.0896 0.9143 0.2579 -0.347 0.728306
beta5 -0.2521 0.7772 0.2605 -0.968 0.333144
cat_pre_wt_log_scale 0.9641 2.6224 0.2556 3.772 0.000162
I(cat_pre_wt_log_scale^2) 0.2450 1.2776 0.1447 1.693 0.090408
var_trtlow_var:cat_pre_wt_log_scale -0.1626 0.8499 0.1998 -0.814 0.415661
beta0:cat_pre_wt_log_scale -0.2230 0.8001 0.2309 -0.966 0.334053
beta5:cat_pre_wt_log_scale -0.4789 0.6195 0.2584 -1.854 0.063807
var_trtlow_var
beta0
beta5
cat_pre_wt_log_scale ***
I(cat_pre_wt_log_scale^2) .
var_trtlow_var:cat_pre_wt_log_scale
beta0:cat_pre_wt_log_scale
beta5:cat_pre_wt_log_scale .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var 0.8945 1.1179 0.5899 1.356
beta0 0.9143 1.0937 0.5515 1.516
beta5 0.7772 1.2867 0.4665 1.295
cat_pre_wt_log_scale 2.6224 0.3813 1.5892 4.328
I(cat_pre_wt_log_scale^2) 1.2776 0.7827 0.9621 1.696
var_trtlow_var:cat_pre_wt_log_scale 0.8499 1.1766 0.5746 1.257
beta0:cat_pre_wt_log_scale 0.8001 1.2499 0.5089 1.258
beta5:cat_pre_wt_log_scale 0.6195 1.6142 0.3734 1.028
Concordance= 0.637 (se = 0.047 )
Likelihood ratio test= 22 on 8 df, p=0.005
Wald test = 22.31 on 8 df, p=0.004
Score (logrank) test = 24.57 on 8 df, p=0.002drop1(surv_m_full, test = "Chisq")Single term deletions
Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + strata(session_id)
Df AIC LRT Pr(>Chi)
<none> 459.50
I(cat_pre_wt_log_scale^2) 1 460.40 2.90 0.08841 .
strata(session_id) 0 787.87 328.37
var_trt:cat_pre_wt_log_scale 1 458.16 0.67 0.41395
beta:cat_pre_wt_log_scale 2 458.98 3.48 0.17532
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1surv_m <- coxph(
Surv(surv_time, observed_dead) ~
(var_trt + beta) + cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
strata(session_id),
data = d %>%
filter(var_trt != "constant") %>%
mutate(
beta = as.factor(as.character(beta))
),
); summary(surv_m)Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) +
cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id),
data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))
n= 128, number of events= 105
coef exp(coef) se(coef) z Pr(>|z|)
var_trtlow_var -0.07077 0.93167 0.20845 -0.340 0.734212
beta0 -0.07410 0.92858 0.25808 -0.287 0.774026
beta5 -0.17722 0.83759 0.25804 -0.687 0.492206
cat_pre_wt_log_scale 0.65947 1.93376 0.18164 3.631 0.000283 ***
I(cat_pre_wt_log_scale^2) 0.29584 1.34425 0.14124 2.095 0.036205 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var 0.9317 1.0733 0.6192 1.402
beta0 0.9286 1.0769 0.5599 1.540
beta5 0.8376 1.1939 0.5051 1.389
cat_pre_wt_log_scale 1.9338 0.5171 1.3545 2.761
I(cat_pre_wt_log_scale^2) 1.3443 0.7439 1.0192 1.773
Concordance= 0.602 (se = 0.047 )
Likelihood ratio test= 18.08 on 5 df, p=0.003
Wald test = 18.29 on 5 df, p=0.003
Score (logrank) test = 19.42 on 5 df, p=0.002drop1(surv_m, test = "Chisq")Single term deletions
Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) + cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + strata(session_id)
Df AIC LRT Pr(>Chi)
<none> 457.42
var_trt 1 455.53 0.12 0.7341150
beta 2 453.89 0.48 0.7883535
cat_pre_wt_log_scale 1 469.05 13.63 0.0002226 ***
I(cat_pre_wt_log_scale^2) 1 459.84 4.42 0.0355376 *
strata(session_id) 0 784.70 327.28
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.4.2 Post-hoc contrasts
2.4.3 Misc stats
r2(surv_m, tolerance = 10^-10) Nagelkerke's R2: 0.1352.4.4 Check model
cox.zph(surv_m) chisq df p
var_trt 1.60 1 0.21
beta 1.43 2 0.49
cat_pre_wt_log_scale 2.61 1 0.11
I(cat_pre_wt_log_scale^2) 1.05 1 0.31
GLOBAL 5.53 5 0.352.4.5 Plot
3 Structural Equation Model
3.1 Sub-models
3.1.1 RGR
subm_rgr <- glmmTMB(
RGR_scale ~
cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq +
sl_mean_obs_log_scale +
prop_explore_logit_scale *
var_toxic_12_scale +
mean_toxic_conc_scale +
area_herb_log_scale +
beta_numeric_scale +
(1|session_id),
data = d2
); summary(subm_rgr) Family: gaussian ( identity )
Formula: RGR_scale ~ cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +
sl_mean_obs_log_scale + prop_explore_logit_scale * var_toxic_12_scale +
mean_toxic_conc_scale + area_herb_log_scale + beta_numeric_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
119.5 149.0 -47.7 95.5 75
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.01245 0.1116
Residual 0.16744 0.4092
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.167
Conditional model:
Estimate Std. Error z value
(Intercept) -0.02766 0.07257 -0.381
cat_pre_wt_log_scale -3.64523 0.41435 -8.797
cat_pre_wt_log_scale_sq -3.01946 0.45628 -6.618
sl_mean_obs_log_scale -0.13590 0.05626 -2.416
prop_explore_logit_scale 0.02370 0.04889 0.485
var_toxic_12_scale -0.01949 0.05115 -0.381
mean_toxic_conc_scale -0.25503 0.07147 -3.568
area_herb_log_scale 0.76522 0.09726 7.868
beta_numeric_scale -0.11265 0.04837 -2.329
prop_explore_logit_scale:var_toxic_12_scale 0.18687 0.04488 4.164
Pr(>|z|)
(Intercept) 0.703071
cat_pre_wt_log_scale < 2e-16 ***
cat_pre_wt_log_scale_sq 3.65e-11 ***
sl_mean_obs_log_scale 0.015710 *
prop_explore_logit_scale 0.627788
var_toxic_12_scale 0.703095
mean_toxic_conc_scale 0.000359 ***
area_herb_log_scale 3.61e-15 ***
beta_numeric_scale 0.019873 *
prop_explore_logit_scale:var_toxic_12_scale 3.13e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.2 Variance in toxin (12 hour)
subm_var_toxic <- glmmTMB(
var_toxic_12_scale ~
sl_mean_obs_log_scale +
beta_numeric_scale +
var_high +
(1|session_id),
data = d2
); summary(subm_var_toxic) Family: gaussian ( identity )
Formula:
var_toxic_12_scale ~ sl_mean_obs_log_scale + beta_numeric_scale +
var_high + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
178.2 193.0 -83.1 166.2 81
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.229e-11 7.892e-06
Residual 3.956e-01 6.290e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.396
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.75075 0.09477 -7.922 2.34e-15 ***
sl_mean_obs_log_scale 0.21267 0.07751 2.744 0.00607 **
beta_numeric_scale -0.15689 0.06811 -2.303 0.02126 *
var_high 1.56021 0.13868 11.250 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.3 Mean step length
subm_sl <- glmmTMB(
sl_mean_obs_log_scale ~
prop_explore_logit_scale +
cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq +
on_toxic_logit_scale +
(1|session_id),
data = d2,
); summary(subm_sl) Family: gaussian ( identity )
Formula:
sl_mean_obs_log_scale ~ prop_explore_logit_scale + cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq + on_toxic_logit_scale + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
223.4 240.7 -104.7 209.4 80
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 4.610e-10 2.147e-05
Residual 6.502e-01 8.064e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.65
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.06163 0.09646 0.639 0.52287
prop_explore_logit_scale -0.26175 0.08987 -2.912 0.00359 **
cat_pre_wt_log_scale 1.37948 0.70121 1.967 0.04915 *
cat_pre_wt_log_scale_sq 1.61095 0.75569 2.132 0.03303 *
on_toxic_logit_scale 0.21435 0.09489 2.259 0.02389 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.4 Mean toxin ingested
subm_toxin_ingested <- glmmTMB(
mean_toxic_conc_scale ~
beta_numeric_scale +
on_toxic_logit_scale +
var_high +
(1|session_id),
data = d2
); summary(subm_toxin_ingested) Family: gaussian ( identity )
Formula:
mean_toxic_conc_scale ~ beta_numeric_scale + on_toxic_logit_scale +
var_high + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
145.7 160.5 -66.8 133.7 81
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.461e-10 2.542e-05
Residual 2.722e-01 5.217e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.272
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.18349 0.07861 2.334 0.0196 *
beta_numeric_scale -0.11683 0.05808 -2.012 0.0443 *
on_toxic_logit_scale 0.28426 0.06393 4.447 8.72e-06 ***
var_high -0.80646 0.11613 -6.945 3.79e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.5 Mean time on toxin
subm_on_toxic <- glmmTMB(
on_toxic_logit_scale ~
var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1|session_id),
family = gaussian(),
data = d2
); summary(subm_on_toxic) Family: gaussian ( identity )
Formula:
on_toxic_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
230.0 247.3 -108.0 216.0 80
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.01667 0.1291
Residual 0.68757 0.8292
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.688
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.22509 0.14140 1.592 0.1114
var_high -0.41809 0.18032 -2.319 0.0204 *
beta_numeric_scale -0.12409 0.09766 -1.271 0.2039
cat_pre_wt_log_scale -0.19306 0.13787 -1.400 0.1614
beta_numeric_scale:cat_pre_wt_log_scale -0.23850 0.11732 -2.033 0.0421 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.6 Area consumption
subm_herb <- glmmTMB(
area_herb_log_scale ~
var_high * cat_pre_wt_log_scale +
(1|session_id),
data = d2,
); summary(subm_herb) Family: gaussian ( identity )
Formula: area_herb_log_scale ~ var_high * cat_pre_wt_log_scale + (1 |
session_id)
Data: d2
AIC BIC logLik deviance df.resid
117.9 132.7 -53.0 105.9 81
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.501e-11 5.917e-06
Residual 1.978e-01 4.448e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.198
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.04339 0.06963 0.623 0.53313
var_high 0.27072 0.10127 2.673 0.00751 **
cat_pre_wt_log_scale 0.55603 0.08378 6.637 3.2e-11 ***
var_high:cat_pre_wt_log_scale -0.26289 0.11558 -2.275 0.02293 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.7 Exploration
subm_exp <- glmmTMB(
prop_explore_logit_scale ~
var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1|session_id),
family = gaussian(),
data = d2
); summary(subm_exp) Family: gaussian ( identity )
Formula:
prop_explore_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
250.1 267.3 -118.0 236.1 80
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.781e-10 1.944e-05
Residual 8.830e-01 9.397e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.883
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.17231 0.14397 -1.197 0.2314
var_high 0.15661 0.20355 0.769 0.4416
beta_numeric_scale -0.14059 0.11040 -1.273 0.2029
cat_pre_wt_log_scale -0.04358 0.12362 -0.353 0.7244
beta_numeric_scale:cat_pre_wt_log_scale 0.25977 0.13242 1.962 0.0498 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2 d-separation test
sem_fit <- psem(
subm_rgr,
subm_herb,
subm_var_toxic,
subm_on_toxic,
subm_toxin_ingested,
subm_exp,
subm_sl,
var_toxic_12_scale %~~% on_toxic_logit_scale, # correlated errors
data = d2
)
sem_summary <- suppressWarnings(suppressMessages(
summary_psem(sem_fit,
no_standardize_x = c(),
.progressBar = FALSE,
direction = "var_toxic_12_scale -> area_herb_log_scale")
))Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.sem_summary$Cstat Fisher.C df P.Value
1 42.469 42 0.451knitr::kable(sem_summary$dTable)| Independ.Claim | Test.Type | DF | Crit.Value | P.Value | |
|---|---|---|---|---|---|
| mean_toxic_conc_scale ~ cat_pre_wt_log_scale + … | coef | 87 | -1.5712 | 0.1161 | |
| var_toxic_12_scale ~ cat_pre_wt_log_scale + … | coef | 87 | -0.1892 | 0.8499 | |
| prop_explore_logit_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | 0.8624 | 0.3884 | |
| mean_toxic_conc_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | 1.8002 | 0.0718 | |
| area_herb_log_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | -0.7057 | 0.4804 | |
| on_toxic_logit_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | 0.7944 | 0.4269 | |
| var_toxic_12_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | 0.3081 | 0.7580 | |
| sl_mean_obs_log_scale ~ beta_numeric_scale + … | coef | 87 | 1.1020 | 0.2704 | |
| area_herb_log_scale ~ beta_numeric_scale + … | coef | 87 | 0.0200 | 0.9840 | |
| RGR_scale ~ var_high + … | coef | 87 | 0.6991 | 0.4845 | |
| sl_mean_obs_log_scale ~ var_high + … | coef | 87 | -1.5883 | 0.1122 | |
| on_toxic_logit_scale ~ RGR_scale + … | coef | 87 | -1.3664 | 0.1718 | |
| mean_toxic_conc_scale ~ sl_mean_obs_log_scale + … | coef | 87 | 0.7459 | 0.4557 | |
| area_herb_log_scale ~ sl_mean_obs_log_scale + … | coef | 87 | 0.8564 | 0.3918 | |
| mean_toxic_conc_scale ~ prop_explore_logit_scale + … | coef | 87 | 1.0591 | 0.2895 | |
| area_herb_log_scale ~ prop_explore_logit_scale + … | coef | 87 | 0.3920 | 0.6951 | |
| on_toxic_logit_scale ~ prop_explore_logit_scale + … | coef | 87 | 0.2361 | 0.8134 | |
| var_toxic_12_scale ~ prop_explore_logit_scale + … | coef | 87 | 1.1110 | 0.2666 | |
| area_herb_log_scale ~ mean_toxic_conc_scale + … | coef | 87 | -0.2457 | 0.8059 | |
| var_toxic_12_scale ~ mean_toxic_conc_scale + … | coef | 87 | -1.3086 | 0.1907 | |
| on_toxic_logit_scale ~ area_herb_log_scale + … | coef | 87 | -0.7213 | 0.4707 |
3.3 Model summary
3.3.1 Model stats
sem_summary$ChiSqNULLsem_summary$AIC AIC AICc K n
1 1264.813 1276.432 51 87sem_summary$R2 Response family link method Marginal Conditional
1 RGR_scale gaussian identity none 0.75 0.77
2 area_herb_log_scale gaussian identity none 0.42 NA
3 var_toxic_12_scale gaussian identity none 0.60 NA
4 on_toxic_logit_scale gaussian identity none 0.19 0.21
5 mean_toxic_conc_scale gaussian identity none 0.55 NA
6 prop_explore_logit_scale gaussian identity none 0.06 NA
7 sl_mean_obs_log_scale gaussian identity none 0.19 NA3.3.2 Standardized coefficients
knitr::kable(
cbind(sem_summary$coefficients,"arrow_size"=(abs(as.numeric(sem_summary$coefficients$Std.Estimate))*20))
)| Response | Predictor | Estimate | Std.Error | DF | Crit.Value | P.Value | Std.Estimate | star | arrow_size |
|---|---|---|---|---|---|---|---|---|---|
| RGR_scale | cat_pre_wt_log_scale | -3.6452 | 0.4144 | 87 | -8.7974 | 0.0000 | -3.5981 | *** | 71.962 |
| RGR_scale | cat_pre_wt_log_scale_sq | -3.0195 | 0.4563 | 87 | -6.6176 | 0.0000 | -2.7708 | *** | 55.416 |
| RGR_scale | sl_mean_obs_log_scale | -0.1359 | 0.0563 | 87 | -2.4156 | 0.0157 | -0.1448 | * | 2.896 |
| RGR_scale | prop_explore_logit_scale | 0.0237 | 0.0489 | 87 | 0.4848 | 0.6278 | 0.0274 | 0.548 | |
| RGR_scale | var_toxic_12_scale | -0.0195 | 0.0511 | 87 | -0.3811 | 0.7031 | -0.0230 | 0.460 | |
| RGR_scale | mean_toxic_conc_scale | -0.2550 | 0.0715 | 87 | -3.5683 | 0.0004 | -0.2353 | *** | 4.706 |
| RGR_scale | area_herb_log_scale | 0.7652 | 0.0973 | 87 | 7.8679 | 0.0000 | 0.5338 | *** | 10.676 |
| RGR_scale | beta_numeric_scale | -0.1126 | 0.0484 | 87 | -2.3287 | 0.0199 | -0.1337 | * | 2.674 |
| RGR_scale | prop_explore_logit_scale:var_toxic_12_scale | 0.1869 | 0.0449 | 87 | 4.1640 | 0.0000 | 0.2280 | *** | 4.560 |
| area_herb_log_scale | var_high | 0.2707 | 0.1013 | 87 | 2.6733 | 0.0075 | 0.2315 | ** | 4.630 |
| area_herb_log_scale | cat_pre_wt_log_scale | 0.5560 | 0.0838 | 87 | 6.6371 | 0.0000 | 0.7868 | *** | 15.736 |
| area_herb_log_scale | var_high:cat_pre_wt_log_scale | -0.2629 | 0.1156 | 87 | -2.2745 | 0.0229 | -0.2795 | * | 5.590 |
| var_toxic_12_scale | sl_mean_obs_log_scale | 0.2127 | 0.0775 | 87 | 2.7437 | 0.0061 | 0.1918 | ** | 3.836 |
| var_toxic_12_scale | beta_numeric_scale | -0.1569 | 0.0681 | 87 | -2.3033 | 0.0213 | -0.1575 | * | 3.150 |
| var_toxic_12_scale | var_high | 1.5602 | 0.1387 | 87 | 11.2503 | 0.0000 | 0.7875 | *** | 15.750 |
| on_toxic_logit_scale | var_high | -0.4181 | 0.1803 | 87 | -2.3186 | 0.0204 | -0.2235 | * | 4.470 |
| on_toxic_logit_scale | beta_numeric_scale | -0.1241 | 0.0977 | 87 | -1.2706 | 0.2039 | -0.1320 | 2.640 | |
| on_toxic_logit_scale | cat_pre_wt_log_scale | -0.1931 | 0.1379 | 87 | -1.4003 | 0.1614 | -0.1708 | 3.416 | |
| on_toxic_logit_scale | beta_numeric_scale:cat_pre_wt_log_scale | -0.2385 | 0.1173 | 87 | -2.0329 | 0.0421 | -0.2116 | * | 4.232 |
| mean_toxic_conc_scale | beta_numeric_scale | -0.1168 | 0.0581 | 87 | -2.0116 | 0.0443 | -0.1503 | * | 3.006 |
| mean_toxic_conc_scale | on_toxic_logit_scale | 0.2843 | 0.0639 | 87 | 4.4467 | 0.0000 | 0.3438 | *** | 6.876 |
| mean_toxic_conc_scale | var_high | -0.8065 | 0.1161 | 87 | -6.9447 | 0.0000 | -0.5214 | *** | 10.428 |
| prop_explore_logit_scale | var_high | 0.1566 | 0.2035 | 87 | 0.7694 | 0.4416 | 0.0808 | 1.616 | |
| prop_explore_logit_scale | beta_numeric_scale | -0.1406 | 0.1104 | 87 | -1.2734 | 0.2029 | -0.1443 | 2.886 | |
| prop_explore_logit_scale | cat_pre_wt_log_scale | -0.0436 | 0.1236 | 87 | -0.3526 | 0.7244 | -0.0372 | 0.744 | |
| prop_explore_logit_scale | beta_numeric_scale:cat_pre_wt_log_scale | 0.2598 | 0.1324 | 87 | 1.9617 | 0.0498 | 0.2224 | * | 4.448 |
| sl_mean_obs_log_scale | prop_explore_logit_scale | -0.2618 | 0.0899 | 87 | -2.9124 | 0.0036 | -0.2840 | ** | 5.680 |
| sl_mean_obs_log_scale | cat_pre_wt_log_scale | 1.3795 | 0.7012 | 87 | 1.9673 | 0.0491 | 1.2778 | * | 25.556 |
| sl_mean_obs_log_scale | cat_pre_wt_log_scale_sq | 1.6109 | 0.7557 | 87 | 2.1318 | 0.0330 | 1.3873 | * | 27.746 |
| sl_mean_obs_log_scale | on_toxic_logit_scale | 0.2144 | 0.0949 | 87 | 2.2589 | 0.0239 | 0.2244 | * | 4.488 |
| ~~var_toxic_12_scale | ~~on_toxic_logit_scale | 0.5323 | - | 87 | 5.7626 | 0.0000 | 0.5323 | *** | 10.646 |
3.3.3 Post hoc contrast
contrast_by_pre_wt(subm_herb, "var_high", "trt")$emmeans
cat_pre_wt_log_scale = -2:
var_high emmean SE df lower.CL upper.CL
0 -1.069 0.200 81 -1.467 -0.670
1 -0.272 0.198 81 -0.666 0.122
cat_pre_wt_log_scale = 2:
var_high emmean SE df lower.CL upper.CL
0 1.155 0.161 81 0.836 1.475
1 0.900 0.150 81 0.603 1.198
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 0.797 0.281 81 2.830 0.0059
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 -0.255 0.219 81 -1.162 0.2485sd(subm_herb$frame$var_high)/sd(subm_herb$frame$area_herb_log_scale)[1] 0.8549509emmeans::emtrends(subm_exp,
~ beta_numeric_scale | cat_pre_wt_log_scale,
var = "beta_numeric_scale",
at = list(cat_pre_wt_log_scale = c(-2,2)))cat_pre_wt_log_scale = -2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.88e-17 -0.660 0.323 80 -1.30 -0.0173
cat_pre_wt_log_scale = 2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.88e-17 0.379 0.246 80 -0.11 0.8677
Results are averaged over the levels of: var_high
Confidence level used: 0.95 sd(subm_exp$frame$beta_numeric_scale) * sd(subm_exp$frame$prop_explore_logit_scale)[1] 0.974435emmeans::emtrends(subm_on_toxic,
~ beta_numeric_scale | cat_pre_wt_log_scale,
var = "beta_numeric_scale",
at = list(cat_pre_wt_log_scale = c(-2,2)))cat_pre_wt_log_scale = -2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.88e-17 0.353 0.286 80 -0.217 0.923
cat_pre_wt_log_scale = 2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.88e-17 -0.601 0.217 80 -1.034 -0.169
Results are averaged over the levels of: var_high
Confidence level used: 0.95 sd(subm_on_toxic$frame$beta_numeric_scale) * sd(subm_on_toxic$frame$on_toxic_logit_scale)[1] 0.94033613.4 SEM plots
knitr::include_graphics(paste0(getwd(), "/graphs/manuscript1_figures/SEM_simplified.png"), rel_path = FALSE)
4 ISSF parameters
4.1 Selection strength state 1
4.1.1 Model summary
m_select_s1_full <- glmmTMB(
s1.less_toxic.est ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_select_s1_full) Family: gaussian ( identity )
Formula:
s1.less_toxic.est ~ cat_pre_wt_log_scale * (beta + var_trt) +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
74.9 98.7 -27.5 54.9 70
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.001837 0.04286
Residual 0.114702 0.33868
Number of obs: 80, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.115
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.0114538 0.0806500 -0.142 0.887
cat_pre_wt_log_scale 0.0836438 0.0902695 0.927 0.354
beta0 0.0852297 0.0934340 0.912 0.362
beta5 0.0170840 0.1040511 0.164 0.870
var_trtlow_var -0.0195540 0.0794866 -0.246 0.806
cat_pre_wt_log_scale:beta0 0.0730532 0.1003755 0.728 0.467
cat_pre_wt_log_scale:beta5 0.0823297 0.1245307 0.661 0.509
cat_pre_wt_log_scale:var_trtlow_var -0.0003578 0.0910877 -0.004 0.997car::Anova(m_select_s1_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: s1.less_toxic.est
Chisq Df Pr(>Chisq)
(Intercept) 0.0202 1 0.8871
cat_pre_wt_log_scale 0.8586 1 0.3541
beta 0.9317 2 0.6276
var_trt 0.0605 1 0.8057
cat_pre_wt_log_scale:beta 0.6707 2 0.7151
cat_pre_wt_log_scale:var_trt 0.0000 1 0.9969m_select_s1 <- glmmTMB(
s1.less_toxic.est ~
cat_pre_wt_log_scale + beta + var_trt +
(1|session_id),
data = d3,
); summary(m_select_s1) Family: gaussian ( identity )
Formula:
s1.less_toxic.est ~ cat_pre_wt_log_scale + beta + var_trt + (1 |
session_id)
Data: d3
AIC BIC logLik deviance df.resid
69.6 86.3 -27.8 55.6 73
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.001818 0.04264
Residual 0.115694 0.34014
Number of obs: 80, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.116
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.01987 0.08008 -0.248 0.80402
cat_pre_wt_log_scale 0.13406 0.04556 2.942 0.00326 **
beta0 0.10143 0.09093 1.115 0.26468
beta5 0.03972 0.09790 0.406 0.68494
var_trtlow_var -0.02615 0.07649 -0.342 0.73245
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_select_s1, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: s1.less_toxic.est
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 8.6583 1 0.003256 **
beta 1.2722 2 0.529358
var_trt 0.1169 1 0.732451
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.1.2 Post-hoc contrasts
contrast_by_pre_wt(m_select_s1, "cat_pre_wt_log_scale")$emmeans
cat_pre_wt_log_scale emmean SE df lower.CL upper.CL
-2 -0.254 0.111 73 -0.475 -0.0331
2 0.282 0.091 73 0.101 0.4636
Results are averaged over the levels of: beta, var_trt
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio
(cat_pre_wt_log_scale-2) - cat_pre_wt_log_scale2 -0.536 0.182 73 -2.943
p.value
0.0044
Results are averaged over the levels of: beta, var_trt 4.1.3 Misc stats
r2(m_select_s1, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.129
Marginal R2: 0.1154.1.4 Check model
check_mod(m_select_s1)
4.1.5 Plot
marginal_effects(m_select_s1, terms = c("cat_pre_wt_log_scale")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(yhat))) +
geom_hline(aes(yintercept = 1), color = "tomato", linetype = "dashed", linewidth = 1) +
geom_ribbon(aes(ymax = exp(upper), ymin = exp(lower)), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_select_s1$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(s1.less_toxic.est)),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10", breaks = c(0.5, 2, 1)) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
labs(x = "Pre-weight (g)", y = "Less toxic diet selection strength (odds ratio)")
4.2 Selection strength state 2
4.2.1 Model summary
m_select_s2_full <- glmmTMB(
s2.less_toxic.est ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_select_s2_full) Family: gaussian ( identity )
Formula:
s2.less_toxic.est ~ cat_pre_wt_log_scale * (beta + var_trt) +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
162.4 186.1 -71.2 142.4 69
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 8.893e-11 9.430e-06
Residual 3.550e-01 5.958e-01
Number of obs: 79, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.355
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.16704 0.13688 1.220 0.222
cat_pre_wt_log_scale 0.10424 0.15818 0.659 0.510
beta0 -0.10199 0.16549 -0.616 0.538
beta5 -0.11331 0.18192 -0.623 0.533
var_trtlow_var -0.02173 0.14187 -0.153 0.878
cat_pre_wt_log_scale:beta0 0.27329 0.18204 1.501 0.133
cat_pre_wt_log_scale:beta5 -0.30413 0.21867 -1.391 0.164
cat_pre_wt_log_scale:var_trtlow_var -0.23822 0.16499 -1.444 0.149# Somehow the wald chisquare test is significant
car::Anova(m_select_s2_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: s2.less_toxic.est
Chisq Df Pr(>Chisq)
(Intercept) 1.4892 1 0.22234
cat_pre_wt_log_scale 0.4343 1 0.50990
beta 0.5224 2 0.77012
var_trt 0.0235 1 0.87825
cat_pre_wt_log_scale:beta 7.5298 2 0.02317 *
cat_pre_wt_log_scale:var_trt 2.0848 1 0.14877
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_select_s2_2 <- glmmTMB(
s2.less_toxic.est ~
cat_pre_wt_log_scale * beta + var_trt +
(1|session_id),
data = d3,
); summary(m_select_s2_2) Family: gaussian ( identity )
Formula:
s2.less_toxic.est ~ cat_pre_wt_log_scale * beta + var_trt + (1 |
session_id)
Data: d3
AIC BIC logLik deviance df.resid
162.4 183.8 -72.2 144.4 70
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 7.426e-11 8.618e-06
Residual 3.643e-01 6.036e-01
Number of obs: 79, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.364
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.18124 0.13832 1.310 0.190
cat_pre_wt_log_scale -0.01899 0.13493 -0.141 0.888
beta0 -0.09318 0.16754 -0.556 0.578
beta5 -0.10415 0.18419 -0.566 0.572
var_trtlow_var -0.08197 0.13738 -0.597 0.551
cat_pre_wt_log_scale:beta0 0.27152 0.18442 1.472 0.141
cat_pre_wt_log_scale:beta5 -0.24385 0.21746 -1.121 0.262# Somehow the wald chisquare test is significant
car::Anova(m_select_s2_2, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: s2.less_toxic.est
Chisq Df Pr(>Chisq)
(Intercept) 1.7170 1 0.19008
cat_pre_wt_log_scale 0.0198 1 0.88809
beta 0.4282 2 0.80726
var_trt 0.3560 1 0.55072
cat_pre_wt_log_scale:beta 6.2016 2 0.04501 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# hmmm actually??
# Drop the interaction and compare with LRT directly
m_select_s2 <- glmmTMB(
s2.less_toxic.est ~
cat_pre_wt_log_scale + var_trt + beta +
(1|session_id),
data = d3,
); summary(m_select_s2) Family: gaussian ( identity )
Formula:
s2.less_toxic.est ~ cat_pre_wt_log_scale + var_trt + beta + (1 |
session_id)
Data: d3
AIC BIC logLik deviance df.resid
164.4 181.0 -75.2 150.4 72
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 1.016e-10 1.008e-05
Residual 3.929e-01 6.269e-01
Number of obs: 79, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.393
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.15638 0.14205 1.101 0.271
cat_pre_wt_log_scale 0.03591 0.08405 0.427 0.669
var_trtlow_var -0.05883 0.14169 -0.415 0.678
beta0 -0.02543 0.16765 -0.152 0.879
beta5 -0.19507 0.17859 -1.092 0.275# LRT
anova(m_select_s2, m_select_s2_2)Data: d3
Models:
m_select_s2: s2.less_toxic.est ~ cat_pre_wt_log_scale + var_trt + beta + (1 | , zi=~0, disp=~1
m_select_s2: session_id), zi=~0, disp=~1
m_select_s2_2: s2.less_toxic.est ~ cat_pre_wt_log_scale * beta + var_trt + (1 | , zi=~0, disp=~1
m_select_s2_2: session_id), zi=~0, disp=~1
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_select_s2 7 164.40 180.99 -75.200 150.40
m_select_s2_2 9 162.43 183.75 -72.215 144.43 5.9702 2 0.05054 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_select_s2, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: s2.less_toxic.est
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 0.1825 1 0.6692
var_trt 0.1724 1 0.6780
beta 1.3823 2 0.50104.2.2 Post-hoc contrasts
emmeans::emmeans(m_select_s2, ~1) 1 emmean SE df lower.CL upper.CL
overall 0.0634 0.0709 72 -0.0779 0.205
Results are averaged over the levels of: var_trt, beta
Confidence level used: 0.95 4.2.3 Misc stats
r2(m_select_s2, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.020
Marginal R2: 0.0204.2.4 Check model
# Diagnostic plot looks better with m_select_s2 than m_select_s2_2
check_mod(m_select_s2)
4.2.5 Plot
marginal_effects(m_select_s2, terms = c("cat_pre_wt_log_scale")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(yhat))) +
geom_hline(aes(yintercept = 1), color = "tomato", linetype = "dashed", linewidth = 1) +
geom_ribbon(aes(ymax = exp(upper), ymin = exp(lower)), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_select_s2$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(s2.less_toxic.est)),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10", breaks = c(0.5, 2, 1)) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
labs(x = "Pre-weight (g)", y = "Less toxic diet selection strength (odds ratio)")
4.3 Arrestment state 1
4.3.1 Model summary
m_arrest_s1_full <- glmmTMB(
log(k1) ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_arrest_s1_full) Family: gaussian ( identity )
Formula:
log(k1) ~ cat_pre_wt_log_scale * (beta + var_trt) + (1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
217.2 241.5 -98.6 197.2 74
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 4.632e-10 2.152e-05
Residual 6.123e-01 7.825e-01
Number of obs: 84, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.612
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.40354 0.17829 2.263 0.0236 *
cat_pre_wt_log_scale 0.12847 0.20516 0.626 0.5312
beta0 -0.32607 0.21241 -1.535 0.1248
beta5 -0.52396 0.23482 -2.231 0.0257 *
var_trtlow_var -0.28762 0.17941 -1.603 0.1089
cat_pre_wt_log_scale:beta0 -0.21810 0.23027 -0.947 0.3436
cat_pre_wt_log_scale:beta5 0.67672 0.28167 2.402 0.0163 *
cat_pre_wt_log_scale:var_trtlow_var 0.07069 0.20685 0.342 0.7325
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_arrest_s1_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(k1)
Chisq Df Pr(>Chisq)
(Intercept) 5.1231 1 0.023609 *
cat_pre_wt_log_scale 0.3921 1 0.531188
beta 5.2345 2 0.073005 .
var_trt 2.5699 1 0.108913
cat_pre_wt_log_scale:beta 11.1326 2 0.003825 **
cat_pre_wt_log_scale:var_trt 0.1168 1 0.732528
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_arrest_s1 <- glmmTMB(
log(k1) ~
cat_pre_wt_log_scale * beta + var_trt +
(1|session_id),
data = d3,
); summary(m_arrest_s1) Family: gaussian ( identity )
Formula:
log(k1) ~ cat_pre_wt_log_scale * beta + var_trt + (1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
215.3 237.2 -98.6 197.3 75
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.684e-10 1.919e-05
Residual 6.132e-01 7.831e-01
Number of obs: 84, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.613
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.3998 0.1781 2.245 0.0248 *
cat_pre_wt_log_scale 0.1651 0.1750 0.944 0.3453
beta0 -0.3331 0.2115 -1.575 0.1153
beta5 -0.5260 0.2349 -2.239 0.0252 *
var_trtlow_var -0.2708 0.1726 -1.568 0.1168
cat_pre_wt_log_scale:beta0 -0.2114 0.2296 -0.921 0.3572
cat_pre_wt_log_scale:beta5 0.6591 0.2771 2.378 0.0174 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_arrest_s1, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(k1)
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 4.6103 1 0.031780 *
beta 3.2439 2 0.197512
var_trt 2.4603 1 0.116755
cat_pre_wt_log_scale:beta 11.3151 2 0.003491 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_arrest_s1, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(k1)
Chisq Df Pr(>Chisq)
(Intercept) 5.0410 1 0.024754 *
cat_pre_wt_log_scale 0.8905 1 0.345338
beta 5.3120 2 0.070229 .
var_trt 2.4603 1 0.116755
cat_pre_wt_log_scale:beta 11.3151 2 0.003491 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.3.2 Post-hoc contrasts
contrast_by_pre_wt(m_arrest_s1, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df lower.CL upper.CL
-5 -0.0658 0.419 75 -0.9000 0.768
0 0.0239 0.349 75 -0.6718 0.719
5 -1.9099 0.535 75 -2.9748 -0.845
cat_pre_wt_log_scale = 2:
beta emmean SE df lower.CL upper.CL
-5 0.5947 0.344 75 -0.0907 1.280
0 -0.1612 0.308 75 -0.7751 0.453
5 1.3869 0.375 75 0.6390 2.135
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.0897 0.546 75 -0.164 0.9852
(beta-5) - beta5 1.8441 0.681 75 2.710 0.0225
beta0 - beta5 1.9338 0.638 75 3.032 0.0093
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 0.7560 0.462 75 1.638 0.2363
(beta-5) - beta5 -0.7922 0.511 75 -1.549 0.2740
beta0 - beta5 -1.5481 0.486 75 -3.186 0.0059
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates emmeans::emtrends(m_arrest_s1, var = "cat_pre_wt_log_scale", specs = "beta") beta cat_pre_wt_log_scale.trend SE df lower.CL upper.CL
-5 0.1651 0.175 75 -0.183 0.514
0 -0.0463 0.148 75 -0.342 0.249
5 0.8242 0.214 75 0.399 1.250
Results are averaged over the levels of: var_trt
Confidence level used: 0.95 4.3.3 Misc stats
r2(m_arrest_s1, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.214
Marginal R2: 0.2144.3.4 Check model
check_mod(m_arrest_s1)
Warning: Maximum value of original data is not included in the
replicated data.
Model may not capture the variation of the data.
4.3.5 Plot
marginal_effects(m_arrest_s1, terms = c("cat_pre_wt_log_scale", "beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat), color = beta)) +
geom_ribbon(aes(ymax = (upper), ymin = lower, fill = beta, color = NULL), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_arrest_s1$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(`log(k1)`)),
size = 3, alpha = 0.5
) +
theme_bw(base_size = 15) +
geom_hline(aes(yintercept = 1), color = "black", linetype = "dashed", linewidth = 1) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("color","fill")) +
labs(x = "Pre-weight (g)", y = "Less toxic diet arresetment strength (ratio)",
color = expression(beta), fill = expression(beta))
4.4 Arrestment state 2
4.4.1 Model summary
m_arrest_s2_full <- glmmTMB(
log(k2) ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_arrest_s2_full) Family: gaussian ( identity )
Formula:
log(k2) ~ cat_pre_wt_log_scale * (beta + var_trt) + (1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
73.7 98.0 -26.9 53.7 74
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.005218 0.07224
Residual 0.107167 0.32736
Number of obs: 84, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.107
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.10940 0.08207 1.333 0.183
cat_pre_wt_log_scale 0.05992 0.09223 0.650 0.516
beta0 -0.07309 0.08965 -0.815 0.415
beta5 0.03636 0.09864 0.369 0.712
var_trtlow_var -0.02911 0.07537 -0.386 0.699
cat_pre_wt_log_scale:beta0 0.08419 0.09647 0.873 0.383
cat_pre_wt_log_scale:beta5 -0.02545 0.11795 -0.216 0.829
cat_pre_wt_log_scale:var_trtlow_var -0.09936 0.08662 -1.147 0.251car::Anova(m_arrest_s2_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(k2)
Chisq Df Pr(>Chisq)
(Intercept) 1.7769 1 0.1825
cat_pre_wt_log_scale 0.4221 1 0.5159
beta 1.4394 2 0.4869
var_trt 0.1492 1 0.6993
cat_pre_wt_log_scale:beta 1.2465 2 0.5362
cat_pre_wt_log_scale:var_trt 1.3157 1 0.2514m_arrest_s2 <- glmmTMB(
log(k2) ~
cat_pre_wt_log_scale + beta + var_trt +
(1|session_id),
data = d3,
); summary(m_arrest_s2) Family: gaussian ( identity )
Formula:
log(k2) ~ cat_pre_wt_log_scale + beta + var_trt + (1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
69.8 86.8 -27.9 55.8 77
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.005022 0.07087
Residual 0.110044 0.33173
Number of obs: 84, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.11
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.10749 0.08197 1.311 0.190
cat_pre_wt_log_scale 0.04245 0.05228 0.812 0.417
beta0 -0.04860 0.08821 -0.551 0.582
beta5 0.03345 0.09193 0.364 0.716
var_trtlow_var -0.05442 0.07288 -0.747 0.455car::Anova(m_arrest_s2, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(k2)
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 0.6592 1 0.4168
beta 0.8555 2 0.6520
var_trt 0.5576 1 0.45524.4.2 Post-hoc contrasts
emmeans::emmeans(m_arrest_s2, ~ 1) 1 emmean SE df lower.CL upper.CL
overall 0.0861 0.0492 77 -0.0118 0.184
Results are averaged over the levels of: beta, var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95 4.4.3 Misc stats
r2(m_arrest_s2, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.073
Marginal R2: 0.0314.4.4 Check model
check_mod(m_arrest_s2)
4.4.5 Plot
marginal_effects(m_arrest_s2, terms = c("cat_pre_wt_log_scale", "beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat), color = beta)) +
geom_ribbon(aes(ymax = (upper), ymin = lower, fill = beta, color = NULL), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_arrest_s2$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(`log(k2)`)),
size = 3, alpha = 0.5
) +
theme_bw(base_size = 15) +
geom_hline(aes(yintercept = 1), color = "black", linetype = "dashed", linewidth = 1) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("color","fill")) +
labs(x = "Pre-weight (g)", y = "Less toxic diet arresetment strength (ratio)",
color = expression(beta), fill = expression(beta))
4.5 Toxic diet step length state 1
4.5.1 Model summary
m_sl_pred_s1_full <- glmmTMB(
log(sl_mean_pred1) ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_sl_pred_s1_full) Family: gaussian ( identity )
Formula:
log(sl_mean_pred1) ~ cat_pre_wt_log_scale * (beta + var_trt) +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
269.0 293.3 -124.5 249.0 74
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.832e-10 2.614e-05
Residual 1.135e+00 1.065e+00
Number of obs: 84, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 1.13
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.22146 0.24273 17.392 <2e-16 ***
cat_pre_wt_log_scale 0.09935 0.27931 0.356 0.722
beta0 -0.12085 0.28919 -0.418 0.676
beta5 0.15813 0.31969 0.495 0.621
var_trtlow_var 0.22314 0.24426 0.914 0.361
cat_pre_wt_log_scale:beta0 0.05220 0.31351 0.167 0.868
cat_pre_wt_log_scale:beta5 -0.27790 0.38348 -0.725 0.469
cat_pre_wt_log_scale:var_trtlow_var 0.08592 0.28162 0.305 0.760
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_sl_pred_s1_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(sl_mean_pred1)
Chisq Df Pr(>Chisq)
(Intercept) 302.4652 1 <2e-16 ***
cat_pre_wt_log_scale 0.1265 1 0.7221
beta 0.8192 2 0.6639
var_trt 0.8345 1 0.3610
cat_pre_wt_log_scale:beta 0.8433 2 0.6560
cat_pre_wt_log_scale:var_trt 0.0931 1 0.7603
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_sl_pred_s1 <- glmmTMB(
log(sl_mean_pred1) ~
cat_pre_wt_log_scale + beta + var_trt +
(1|session_id),
data = d3,
); summary(m_sl_pred_s1) Family: gaussian ( identity )
Formula: log(sl_mean_pred1) ~ cat_pre_wt_log_scale + beta + var_trt +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
264.2 281.2 -125.1 250.2 77
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 5.172e-10 2.274e-05
Residual 1.151e+00 1.073e+00
Number of obs: 84, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 1.15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.21702 0.24165 17.451 <2e-16 ***
cat_pre_wt_log_scale 0.10543 0.13672 0.771 0.441
beta0 -0.12393 0.28269 -0.438 0.661
beta5 0.04361 0.29592 0.147 0.883
var_trtlow_var 0.26089 0.23477 1.111 0.266
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_sl_pred_s1, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(sl_mean_pred1)
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 0.5947 1 0.4406
beta 0.3749 2 0.8291
var_trt 1.2349 1 0.26654.5.2 Post-hoc contrasts
emmeans::emmeans(m_sl_pred_s1, ~ 1) 1 emmean SE df lower.CL upper.CL
overall 4.35 0.117 77 4.11 4.58
Results are averaged over the levels of: beta, var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95 4.5.3 Misc stats
r2(m_sl_pred_s1, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.026
Marginal R2: 0.0264.5.4 Check model
check_mod(m_sl_pred_s1)
4.5.5 Plot
marginal_effects(m_sl_pred_s1, terms = c("cat_pre_wt_log_scale")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat) / 1000 * 12)) +
geom_ribbon(aes(ymax = (upper)/ 1000 * 12, ymin = lower/ 1000 * 12), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_sl_pred_s1$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale),
y = exp(`log(sl_mean_pred1)`) / 1000 * 12),
size = 3, alpha = 0.5
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10", labels = fancy_scientific) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("color","fill")) +
labs(x = "Cat pre-weight (g)", y = "Selection free step length \non toxic diet (cm)")
4.6 Toxic diet step length state 2
4.6.1 Model summary
m_sl_pred_s2_full <- glmmTMB(
log(sl_mean_pred3) ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_sl_pred_s2_full) Family: gaussian ( identity )
Formula:
log(sl_mean_pred3) ~ cat_pre_wt_log_scale * (beta + var_trt) +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
56.1 80.4 -18.0 36.1 74
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.709e-12 2.59e-06
Residual 8.997e-02 3.00e-01
Number of obs: 84, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.09
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.92286 0.06834 28.136 <2e-16 ***
cat_pre_wt_log_scale 0.05769 0.07864 0.734 0.463
beta0 -0.05293 0.08142 -0.650 0.516
beta5 0.02926 0.09001 0.325 0.745
var_trtlow_var -0.02306 0.06877 -0.335 0.737
cat_pre_wt_log_scale:beta0 0.05333 0.08827 0.604 0.546
cat_pre_wt_log_scale:beta5 -0.11443 0.10797 -1.060 0.289
cat_pre_wt_log_scale:var_trtlow_var 0.01622 0.07929 0.205 0.838
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_sl_pred_s2_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(sl_mean_pred3)
Chisq Df Pr(>Chisq)
(Intercept) 791.6069 1 <2e-16 ***
cat_pre_wt_log_scale 0.5382 1 0.4632
beta 0.9686 2 0.6161
var_trt 0.1124 1 0.7374
cat_pre_wt_log_scale:beta 2.6393 2 0.2672
cat_pre_wt_log_scale:var_trt 0.0419 1 0.8379
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_sl_pred_s2 <- glmmTMB(
log(sl_mean_pred3) ~
cat_pre_wt_log_scale + beta + var_trt +
(1|session_id),
data = d3,
); summary(m_sl_pred_s2) Family: gaussian ( identity )
Formula: log(sl_mean_pred3) ~ cat_pre_wt_log_scale + beta + var_trt +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
53.1 70.1 -19.6 39.1 77
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 5.339e-12 2.311e-06
Residual 9.327e-02 3.054e-01
Number of obs: 84, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.0933
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.91914 0.06880 27.896 <2e-16 ***
cat_pre_wt_log_scale 0.06490 0.03892 1.667 0.0955 .
beta0 -0.04648 0.08048 -0.577 0.5636
beta5 -0.01791 0.08425 -0.213 0.8316
var_trtlow_var -0.01312 0.06684 -0.196 0.8443
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_sl_pred_s2, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(sl_mean_pred3)
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 2.7798 1 0.09546 .
beta 0.3419 2 0.84285
var_trt 0.0385 1 0.84435
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.6.2 Post-hoc contrasts
emmeans::emmeans(m_sl_pred_s2, ~ 1) 1 emmean SE df lower.CL upper.CL
overall 1.91 0.0334 77 1.84 1.97
Results are averaged over the levels of: beta, var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95 4.6.3 Misc stats
r2(m_sl_pred_s2, tolerance = 10^-10)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.0384.6.4 Check model
check_mod(m_sl_pred_s2)
4.6.5 Plot
marginal_effects(m_sl_pred_s2, terms = c("cat_pre_wt_log_scale")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat) / 1000 * 12)) +
geom_ribbon(aes(ymax = (upper)/ 1000 * 12, ymin = lower/ 1000 * 12), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_sl_pred_s2$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale),
y = exp(`log(sl_mean_pred3)`) / 1000 * 12),
size = 3, alpha = 0.5
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("color","fill")) +
labs(x = "Pre-weight (g)", y = "Selection free step length \non toxic diet (cm)")
5 ISSF simulation
5.1 Group means
out <- read_csv("simulation/move_rules_sim.csv")
out <- out %>%
mutate(
scale = round(scale / 1000 * 12, 2),
rss = as.factor(round(exp(rss), 2))
)
out %>%
group_by(b, rss, scale, k) %>%
summarise(
mean = mean(on_toxic), se = se(on_toxic), sd = sd(on_toxic),
lower = quantile(on_toxic, probs = 0.025), upper = quantile(on_toxic, probs = 0.975)
) %>%
mutate(
z = (mean - 0.5) / se,
P = pnorm(abs(z), lower.tail = FALSE) * 2
) %>%
knitr::kable()| b | rss | scale | k | mean | se | sd | lower | upper | z | P |
|---|---|---|---|---|---|---|---|---|---|---|
| -5 | 0.75 | 0.1 | 0.2 | 0.9424251 | 0.0031393 | 0.0627853 | 0.7781469 | 0.9990010 | 140.9326886 | 0.0000000 |
| -5 | 0.75 | 0.1 | 1.0 | 0.6173901 | 0.0067163 | 0.1343258 | 0.3436563 | 0.8482517 | 17.4784131 | 0.0000000 |
| -5 | 0.75 | 0.1 | 5.0 | 0.1492183 | 0.0026313 | 0.0526266 | 0.0588911 | 0.2557443 | -133.3096814 | 0.0000000 |
| -5 | 0.75 | 0.3 | 0.2 | 0.9423227 | 0.0015801 | 0.0316014 | 0.8631369 | 0.9930070 | 279.9387133 | 0.0000000 |
| -5 | 0.75 | 0.3 | 1.0 | 0.6284740 | 0.0036121 | 0.0722421 | 0.4764236 | 0.7653846 | 35.5676287 | 0.0000000 |
| -5 | 0.75 | 0.3 | 5.0 | 0.3228646 | 0.0020669 | 0.0413379 | 0.2416583 | 0.4056194 | -85.7011511 | 0.0000000 |
| -5 | 0.75 | 1.0 | 0.2 | 0.8792682 | 0.0015648 | 0.0312959 | 0.8191309 | 0.9410589 | 242.3755402 | 0.0000000 |
| -5 | 0.75 | 1.0 | 1.0 | 0.5974151 | 0.0016296 | 0.0325923 | 0.5333666 | 0.6564186 | 59.7779834 | 0.0000000 |
| -5 | 0.75 | 1.0 | 5.0 | 0.4907093 | 0.0013177 | 0.0263531 | 0.4395604 | 0.5394855 | -7.0509466 | 0.0000000 |
| -5 | 1 | 0.1 | 0.2 | 0.9122977 | 0.0047118 | 0.0942358 | 0.6592158 | 1.0000000 | 87.5034337 | 0.0000000 |
| -5 | 1 | 0.1 | 1.0 | 0.4981344 | 0.0069807 | 0.1396140 | 0.2417333 | 0.7634366 | -0.2672560 | 0.7892720 |
| -5 | 1 | 0.1 | 5.0 | 0.1072278 | 0.0020522 | 0.0410437 | 0.0368881 | 0.1968032 | -191.3921273 | 0.0000000 |
| -5 | 1 | 0.3 | 0.2 | 0.9228122 | 0.0020241 | 0.0404811 | 0.8271728 | 0.9830170 | 208.8938658 | 0.0000000 |
| -5 | 1 | 0.3 | 1.0 | 0.4987537 | 0.0037972 | 0.0759442 | 0.3552448 | 0.6473526 | -0.3282027 | 0.7427584 |
| -5 | 1 | 0.3 | 5.0 | 0.2349500 | 0.0017563 | 0.0351252 | 0.1668332 | 0.3037213 | -150.9173885 | 0.0000000 |
| -5 | 1 | 1.0 | 0.2 | 0.8255120 | 0.0020309 | 0.0406185 | 0.7472278 | 0.9111139 | 160.2775744 | 0.0000000 |
| -5 | 1 | 1.0 | 1.0 | 0.4997777 | 0.0016351 | 0.0327023 | 0.4374875 | 0.5654595 | -0.1359402 | 0.8918686 |
| -5 | 1 | 1.0 | 5.0 | 0.3978771 | 0.0012370 | 0.0247397 | 0.3505744 | 0.4475524 | -82.5578881 | 0.0000000 |
| -5 | 1.33 | 0.1 | 0.2 | 0.8868731 | 0.0056609 | 0.1132184 | 0.5920080 | 0.9990260 | 68.3410293 | 0.0000000 |
| -5 | 1.33 | 0.1 | 1.0 | 0.3757318 | 0.0069049 | 0.1380971 | 0.0898851 | 0.6434316 | -17.9972213 | 0.0000000 |
| -5 | 1.33 | 0.1 | 5.0 | 0.0739061 | 0.0014982 | 0.0299632 | 0.0289461 | 0.1390110 | -284.4116825 | 0.0000000 |
| -5 | 1.33 | 0.3 | 0.2 | 0.8830819 | 0.0029752 | 0.0595030 | 0.7412338 | 0.9700549 | 128.7605166 | 0.0000000 |
| -5 | 1.33 | 0.3 | 1.0 | 0.3762687 | 0.0035496 | 0.0709928 | 0.2507243 | 0.5266234 | -34.8573909 | 0.0000000 |
| -5 | 1.33 | 0.3 | 5.0 | 0.1675774 | 0.0014975 | 0.0299509 | 0.1138861 | 0.2278222 | -221.9782717 | 0.0000000 |
| -5 | 1.33 | 1.0 | 0.2 | 0.7501798 | 0.0025080 | 0.0501600 | 0.6542957 | 0.8472028 | 99.7526681 | 0.0000000 |
| -5 | 1.33 | 1.0 | 1.0 | 0.3998776 | 0.0016647 | 0.0332946 | 0.3375874 | 0.4605644 | -60.1433455 | 0.0000000 |
| -5 | 1.33 | 1.0 | 5.0 | 0.3147153 | 0.0010988 | 0.0219754 | 0.2736763 | 0.3556943 | -168.6293330 | 0.0000000 |
| 0 | 0.75 | 0.1 | 0.2 | 0.9337862 | 0.0044369 | 0.0887383 | 0.7415834 | 1.0000000 | 97.7675007 | 0.0000000 |
| 0 | 0.75 | 0.1 | 1.0 | 0.6214011 | 0.0085684 | 0.1713690 | 0.2866384 | 0.9092657 | 14.1683889 | 0.0000000 |
| 0 | 0.75 | 0.1 | 5.0 | 0.1390584 | 0.0032286 | 0.0645729 | 0.0369131 | 0.2767483 | -111.7934353 | 0.0000000 |
| 0 | 0.75 | 0.3 | 0.2 | 0.9485090 | 0.0017937 | 0.0358748 | 0.8560689 | 0.9950050 | 250.0411941 | 0.0000000 |
| 0 | 0.75 | 0.3 | 1.0 | 0.6268906 | 0.0047512 | 0.0950248 | 0.4184316 | 0.7812687 | 26.7068361 | 0.0000000 |
| 0 | 0.75 | 0.3 | 5.0 | 0.2686239 | 0.0025234 | 0.0504690 | 0.1625874 | 0.3776474 | -91.6903984 | 0.0000000 |
| 0 | 0.75 | 1.0 | 0.2 | 0.8988037 | 0.0016735 | 0.0334706 | 0.8261489 | 0.9580420 | 238.3006229 | 0.0000000 |
| 0 | 0.75 | 1.0 | 1.0 | 0.6091533 | 0.0021183 | 0.0423660 | 0.5173826 | 0.6903347 | 51.5287800 | 0.0000000 |
| 0 | 0.75 | 1.0 | 5.0 | 0.4610764 | 0.0015692 | 0.0313849 | 0.3965784 | 0.5195305 | -24.8040389 | 0.0000000 |
| 0 | 1 | 0.1 | 0.2 | 0.9022627 | 0.0056081 | 0.1121629 | 0.6373876 | 1.0000000 | 71.7282847 | 0.0000000 |
| 0 | 1 | 0.1 | 1.0 | 0.5022977 | 0.0085353 | 0.1707059 | 0.1944555 | 0.8171828 | 0.2692001 | 0.7877757 |
| 0 | 1 | 0.1 | 5.0 | 0.0890809 | 0.0022235 | 0.0444710 | 0.0239261 | 0.1888611 | -184.8032596 | 0.0000000 |
| 0 | 1 | 0.3 | 0.2 | 0.9232068 | 0.0025534 | 0.0510686 | 0.7851149 | 0.9950300 | 165.7404909 | 0.0000000 |
| 0 | 1 | 0.3 | 1.0 | 0.4948027 | 0.0047967 | 0.0959345 | 0.2896104 | 0.6853147 | -1.0835101 | 0.2785821 |
| 0 | 1 | 0.3 | 5.0 | 0.1968007 | 0.0020092 | 0.0401848 | 0.1228272 | 0.2748002 | -150.9023306 | 0.0000000 |
| 0 | 1 | 1.0 | 0.2 | 0.8538861 | 0.0022148 | 0.0442958 | 0.7671329 | 0.9390859 | 159.7830069 | 0.0000000 |
| 0 | 1 | 1.0 | 1.0 | 0.5026124 | 0.0021241 | 0.0424812 | 0.4275475 | 0.5865385 | 1.2299033 | 0.2187333 |
| 0 | 1 | 1.0 | 5.0 | 0.3668606 | 0.0015250 | 0.0305002 | 0.3025475 | 0.4245754 | -87.3040484 | 0.0000000 |
| 0 | 1.33 | 0.1 | 0.2 | 0.8764086 | 0.0068052 | 0.1361033 | 0.4984765 | 1.0000000 | 55.3121967 | 0.0000000 |
| 0 | 1.33 | 0.1 | 1.0 | 0.3752597 | 0.0084105 | 0.1682091 | 0.1037962 | 0.7582418 | -14.8315668 | 0.0000000 |
| 0 | 1.33 | 0.1 | 5.0 | 0.0640634 | 0.0016992 | 0.0339841 | 0.0139860 | 0.1429820 | -256.5533298 | 0.0000000 |
| 0 | 1.33 | 0.3 | 0.2 | 0.8817607 | 0.0036300 | 0.0726000 | 0.7192058 | 0.9930070 | 105.1681861 | 0.0000000 |
| 0 | 1.33 | 0.3 | 1.0 | 0.3808516 | 0.0045535 | 0.0910693 | 0.2167333 | 0.5704296 | -26.1665382 | 0.0000000 |
| 0 | 1.33 | 0.3 | 5.0 | 0.1357892 | 0.0016344 | 0.0326884 | 0.0718781 | 0.1988012 | -222.8377936 | 0.0000000 |
| 0 | 1.33 | 1.0 | 0.2 | 0.7972378 | 0.0029964 | 0.0599279 | 0.6683317 | 0.9161339 | 99.1985180 | 0.0000000 |
| 0 | 1.33 | 1.0 | 1.0 | 0.3928521 | 0.0022280 | 0.0445600 | 0.3156094 | 0.4785215 | -48.0915305 | 0.0000000 |
| 0 | 1.33 | 1.0 | 5.0 | 0.2821778 | 0.0013950 | 0.0278996 | 0.2187812 | 0.3337163 | -156.1473668 | 0.0000000 |
| 5 | 0.75 | 0.1 | 0.2 | 0.8456169 | 0.0119521 | 0.2390421 | 0.0009990 | 1.0000000 | 28.9168206 | 0.0000000 |
| 5 | 0.75 | 0.1 | 1.0 | 0.6022253 | 0.0143786 | 0.2875718 | 0.0009990 | 1.0000000 | 7.1095478 | 0.0000000 |
| 5 | 0.75 | 0.1 | 5.0 | 0.1246553 | 0.0055669 | 0.1113376 | 0.0039710 | 0.4073676 | -67.4245913 | 0.0000000 |
| 5 | 0.75 | 0.3 | 0.2 | 0.9369056 | 0.0041112 | 0.0822243 | 0.7109141 | 1.0000000 | 106.2716509 | 0.0000000 |
| 5 | 0.75 | 0.3 | 1.0 | 0.6151499 | 0.0094377 | 0.1887542 | 0.1916833 | 0.9100899 | 12.2010383 | 0.0000000 |
| 5 | 0.75 | 0.3 | 5.0 | 0.1736838 | 0.0037518 | 0.0750362 | 0.0477273 | 0.3277223 | -86.9757179 | 0.0000000 |
| 5 | 0.75 | 1.0 | 0.2 | 0.9375774 | 0.0020378 | 0.0407562 | 0.8341409 | 0.9920330 | 214.7293117 | 0.0000000 |
| 5 | 0.75 | 1.0 | 1.0 | 0.6167907 | 0.0039535 | 0.0790692 | 0.4505245 | 0.7603397 | 29.5413937 | 0.0000000 |
| 5 | 0.75 | 1.0 | 5.0 | 0.3530644 | 0.0032378 | 0.0647564 | 0.2187812 | 0.4565934 | -45.3810152 | 0.0000000 |
| 5 | 1 | 0.1 | 0.2 | 0.8322652 | 0.0120387 | 0.2407742 | 0.0840659 | 1.0000000 | 27.5997408 | 0.0000000 |
| 5 | 1 | 0.1 | 1.0 | 0.5060365 | 0.0147005 | 0.2940094 | 0.0029720 | 0.9980270 | 0.4106307 | 0.6813433 |
| 5 | 1 | 0.1 | 5.0 | 0.0891633 | 0.0041138 | 0.0822755 | 0.0029970 | 0.3087912 | -99.8685114 | 0.0000000 |
| 5 | 1 | 0.3 | 0.2 | 0.9041209 | 0.0054638 | 0.1092765 | 0.6301948 | 1.0000000 | 73.9630157 | 0.0000000 |
| 5 | 1 | 0.3 | 1.0 | 0.5129121 | 0.0096690 | 0.1933807 | 0.1648102 | 0.8911588 | 1.3354058 | 0.1817436 |
| 5 | 1 | 0.3 | 5.0 | 0.1183641 | 0.0026361 | 0.0527216 | 0.0259740 | 0.2200300 | -144.7739158 | 0.0000000 |
| 5 | 1 | 1.0 | 0.2 | 0.9058317 | 0.0028938 | 0.0578754 | 0.7602148 | 0.9880370 | 140.2432388 | 0.0000000 |
| 5 | 1 | 1.0 | 1.0 | 0.5012537 | 0.0043723 | 0.0874459 | 0.3315435 | 0.6896354 | 0.2867480 | 0.7743053 |
| 5 | 1 | 1.0 | 5.0 | 0.2577822 | 0.0027737 | 0.0554733 | 0.1417083 | 0.3526973 | -87.3277366 | 0.0000000 |
| 5 | 1.33 | 0.1 | 0.2 | 0.8062413 | 0.0128033 | 0.2560660 | 0.0899101 | 1.0000000 | 23.9189306 | 0.0000000 |
| 5 | 1.33 | 0.1 | 1.0 | 0.4136064 | 0.0137456 | 0.2749119 | 0.0109391 | 1.0000000 | -6.2851854 | 0.0000000 |
| 5 | 1.33 | 0.1 | 5.0 | 0.0584565 | 0.0028189 | 0.0563777 | 0.0009990 | 0.1988761 | -156.6377091 | 0.0000000 |
| 5 | 1.33 | 0.3 | 0.2 | 0.8744181 | 0.0065334 | 0.1306674 | 0.5342657 | 1.0000000 | 57.3085721 | 0.0000000 |
| 5 | 1.33 | 0.3 | 1.0 | 0.3822128 | 0.0086708 | 0.1734157 | 0.0898851 | 0.7456294 | -13.5843809 | 0.0000000 |
| 5 | 1.33 | 0.3 | 5.0 | 0.0760165 | 0.0019876 | 0.0397526 | 0.0159341 | 0.1638861 | -213.3113441 | 0.0000000 |
| 5 | 1.33 | 1.0 | 0.2 | 0.8645804 | 0.0040493 | 0.0809870 | 0.6832917 | 0.9800699 | 90.0343512 | 0.0000000 |
| 5 | 1.33 | 1.0 | 1.0 | 0.3792957 | 0.0041334 | 0.0826670 | 0.2037712 | 0.5447053 | -29.2025221 | 0.0000000 |
| 5 | 1.33 | 1.0 | 5.0 | 0.1913511 | 0.0024096 | 0.0481928 | 0.0908591 | 0.2747752 | -128.0892573 | 0.0000000 |
5.2 Some fun
# Precompute some summarize grid
z <- out %>%
filter(b != 0) %>%
filter(scale == 1) %>%
group_by(b, rss, k) %>%
dplyr::select(on_toxic) %>%
mutate(
rss = as.numeric(as.character(rss)),
id = seq_along(on_toxic)
)
# Some simple wrapper for hypothesis testing later in the model scenarios
compute_delta <- function(hypothesis, res){
hypothesis %>%
left_join(res, by = c("k", "rss","b"), relationship = "many-to-many") %>%
group_by(cat_size, id) %>%
summarise(
delta = diff(on_toxic)
) %>%
group_by(cat_size) %>%
do(as.data.frame(t(summarise_vec(.$delta)))) %>%
arrange(rev(cat_size))
}
# Get posterior draws
epred_draws <- function(model, newdata){
res <- brms::posterior_epred(model,
newdata = newdata,
allow_new_levels = TRUE,
re_formula = NA) %>%
t()
colnames(res) <- paste0("draw_", seq_len(ncol(res)))
cbind(newdata,res) %>%
gather(key = draw, value = val, contains("draw_"))
}5.3 Model scenarios
# No arrestment, no selection
data.frame("k" = c(1, 1, 1, 1),
"rss" = c(1,1, 1, 1),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z) # A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.00148 0.00150 0.0936 -0.170 0.185
2 l 0.00148 0.00150 0.0936 -0.170 0.185# Size dependent arrestment, no selection
data.frame("k" = c(0.2, 0.2, 5, 5),
"rss" = c(1,1, 1, 1),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z) # A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.0803 0.0889 0.0709 -0.0660 0.199
2 l -0.140 -0.137 0.0623 -0.266 -0.0299# No arrestment, size dependent selection
data.frame("k" = c(1, 1, 1, 1),
"rss" = c(0.75,0.75, 1.33, 1.33),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z)# A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.0194 0.0255 0.0863 -0.160 0.179
2 l -0.0206 -0.0225 0.0891 -0.201 0.173# Size dependent arrestment, size dependent selection
data.frame("k" = c(0.2, 0.2, 5, 5),
"rss" = c(0.75,0.75, 1.33, 1.33),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z) # A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.0583 0.0639 0.0518 -0.0482 0.147
2 l -0.123 -0.121 0.0537 -0.231 -0.0310# Size dependent arrestment in clustered treatment, no selection
data.frame("k" = c(1, 0.2, 1, 5),
"rss" = c(1,1, 1, 1),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z) # A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.406 0.414 0.0639 0.273 0.509
2 l -0.242 -0.239 0.0647 -0.382 -0.127# Size dependent arrestment in clustered treatment, size dependent selection
data.frame("k" = c(1, 0.2, 1, 5),
"rss" = c(0.75,0.75, 1.33, 1.33),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z) # A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.340 0.344 0.0517 0.228 0.428
2 l -0.209 -0.203 0.0588 -0.330 -0.09395.4 Observed means
5.4.1 Summary
# m <- brm(
# on_toxic ~
# beta * cat_pre_wt_log_scale + var_trt +
# (1|session_id),
# family = Beta(),
# data = d %>%
# filter(var_trt != "constant"),
# chains = 4,
# cores = 4,
# prior = c(
# set_prior("normal(0,1.5)", class = "b"),
# set_prior("normal(0,1)", class = "Intercept"),
# set_prior("cauchy(0,0.5)", class = "sd"),
# set_prior("gamma(0.01, 0.01)", class = "phi")
# ),
# iter = 4000,
# control = list(adapt_delta = 0.95)
# )
# saveRDS(m, "invisible/fitted_models/on_toxic_brm.rds")
m <- readRDS("invisible/fitted_models/on_toxic_brm.rds")
m Family: beta
Links: mu = logit; phi = identity
Formula: on_toxic ~ beta * cat_pre_wt_log_scale + var_trt + (1 | session_id)
Data: d %>% filter(var_trt != "constant") (Number of observations: 97)
Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
total post-warmup draws = 8000
Group-Level Effects:
~session_id (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.24 0.19 0.01 0.72 1.00 2231 3150
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept -0.73 0.23 -1.19 -0.27 1.00 4467
beta0 0.14 0.22 -0.29 0.56 1.00 7373
beta5 -0.20 0.24 -0.66 0.28 1.00 6607
cat_pre_wt_log_scale 0.10 0.20 -0.28 0.49 1.00 4581
var_trtlow_var 0.36 0.18 0.01 0.72 1.00 8461
beta0:cat_pre_wt_log_scale -0.30 0.24 -0.76 0.18 1.00 5456
beta5:cat_pre_wt_log_scale -0.72 0.29 -1.29 -0.16 1.00 5511
Tail_ESS
Intercept 4001
beta0 5574
beta5 5527
cat_pre_wt_log_scale 5668
var_trtlow_var 5568
beta0:cat_pre_wt_log_scale 5516
beta5:cat_pre_wt_log_scale 6012
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 4.72 0.66 3.51 6.09 1.00 6699 5717
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).5.4.2 Check model
brms:::plot.brmsfit(m)
brms::pp_check(m,ndraws = 100)
5.4.3 Plot
brms::conditional_effects(m,effects = c("cat_pre_wt_log_scale:beta"), re_formula = NA)[[1]] %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = estimate__, color = beta)) +
geom_ribbon(aes(ymax = upper__, ymin = lower__, color = NULL, fill = beta), alpha = 0.2) +
geom_line(size = 2) +
geom_point(
data = m$data,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = on_toxic),
size = 5,
alpha = 0.5
) +
theme_bw(base_size = 15) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("fill", "color")) +
labs(x = "Pre-weight (g)", y = "Proprtion of time on toxic diet", color = expression(beta), fill = expression(beta))
5.4.4 Observed effect sizes
# Compute observed effect size
epred_draws(m, expand.grid("cat_pre_wt_log_scale" = c(-2, 2),
"beta" = c(-5, 5),
"var_trt" = c("low_var","high_var"),
"session_id" = "foo")) %>%
group_by(cat_pre_wt_log_scale, beta, draw) %>%
summarise(val = mean(val)) %>%
group_by(
cat_pre_wt_log_scale, draw
) %>%
summarise(
delta = diff(val)
) %>%
group_by(cat_pre_wt_log_scale) %>%
do(as.data.frame(t(summarise_vec(.$delta))))# A tibble: 2 × 6
# Groups: cat_pre_wt_log_scale [2]
cat_pre_wt_log_scale mean median sd lower upper
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 -2 0.280 0.284 0.151 -0.0251 0.563
2 2 -0.285 -0.284 0.0970 -0.479 -0.09915.5 Plots
5.5.1 Simplified plot
out <- read_csv("simulation/move_rules_sim.csv")
out <- out %>%
mutate(
scale = round(scale / 1000 * 12, 2),
rss = as.factor(round(exp(rss), 2))
)
out %>%
filter(b != 0) %>%
filter(scale == 1) %>%
group_by(b, rss, k) %>%
dplyr::select(on_toxic) %>%
mutate(
id = seq_along(on_toxic)
) %>%
arrange(b, rss, k, id) %>%
group_by(rss, k, id) %>%
summarise(
delta = diff(on_toxic)
) %>%
ggplot(aes(x = as.factor(k), y = delta, color = factor(rss))) +
geom_hline(aes(yintercept = 0), color = "grey",
size = 1, linetype = "dashed") +
geom_point(position = position_jitterdodge(jitter.width = 0.2,
jitter.height = 0,
dodge.width = 0.5),
alpha = 0.2,
size = 2) +
geom_point(stat = "summary",
position = position_dodge(width = 0.5),
shape = "—",
size = 4,
color = "black",
aes(group = factor(rss))) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
guides(colour = guide_legend(
override.aes = list(alpha = 1, size = 4),
title.position = "top"
)) +
labs(x = "Less toxic diet arresetment strength (ratio)",
y = expression(atop(Change~"in"~proprtion~time~on, paste("more toxic diet",~(beta[5]-beta[-5])))),
color = "Less toxic diet selection strength (odds)") +
scale_color_brewer(type = "seq") -> g;g
5.5.2 Full result
out <- read_csv("simulation/move_rules_sim.csv")
out <- out %>%
mutate(
scale = round(scale / 1000 * 12, 2),
rss = as.factor(round(exp(rss), 2))
)
k_lab <- unique(out$k)
names(k_lab) = sprintf("k = %s", k_lab)
scale_lab <- unique(out$scale)
names(scale_lab) = sprintf("scale = %s", scale_lab)
rss_lab <- unique(out$rss)
names(rss_lab) = sprintf("rss = %s", rss_lab)
out %>%
ggplot(aes(x = as.factor(k), y = on_toxic, color = factor(b))) +
geom_point(position = position_jitterdodge(jitter.width = 0.2,
jitter.height = 0,
dodge.width = 0.5),
alpha = 0.2,
size = 1) +
geom_point(stat = "summary",
position = position_dodge(width = 0.5),
shape = "—",
size = 3,
color = "black",
aes(group = factor(b))) +
theme_bw(base_size = 15) +
guides(colour = guide_legend(
override.aes = list(alpha = 1, size = 3),
)) +
facet_grid(scale ~ rss, labeller = labeller(
rss = reverse_names(rss_lab), scale = reverse_names(scale_lab)
)) +
guides(colour = guide_legend(override.aes = list(alpha = 1))) +
labs(x = "Less toxic diet arresetment strength (ratio)",
y = "Proportion time on toxic diet",
color = expression(beta)) +
scale_color_brewer(type = "qual")->g;g
6 Session Info
sessionInfo()R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8
[2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] spat1f_0.0.1 herbivar_0.2.1 forcats_0.5.1 stringr_1.5.0
[5] dplyr_1.1.2 purrr_1.0.1 readr_2.1.2 tidyr_1.3.0
[9] tibble_3.2.1 tidyverse_1.3.1 extraDistr_1.9.1 imager_0.42.13
[13] magrittr_2.0.3 doSNOW_1.0.20 snow_0.4-4 iterators_1.0.14
[17] foreach_1.5.2 glmmTMB_1.1.7 tidybayes_3.0.2 performance_0.10.4
[21] sjPlot_2.8.15 ggpubr_0.6.0 ggplot2_3.4.2 piecewiseSEM_2.3.0
[25] survival_3.3-1
loaded via a namespace (and not attached):
[1] fs_1.5.2 matrixStats_1.2.0 bitops_1.0-7
[4] sf_1.0-8 devtools_2.4.5 lubridate_1.8.0
[7] httr_1.4.3 RColorBrewer_1.1-3 insight_0.19.3
[10] doParallel_1.0.17 numDeriv_2016.8-1.1 profvis_0.3.7
[13] tools_4.3.1 backports_1.4.1 DT_0.23
[16] sjlabelled_1.2.0 utf8_1.2.2 R6_2.5.1
[19] mgcv_1.8-40 ggdist_3.1.1 urlchecker_1.0.1
[22] withr_2.5.0 sp_1.5-0 readbitmap_0.1.5
[25] gridExtra_2.3 Brobdingnag_1.2-9 prettyunits_1.1.1
[28] cli_3.6.1 shinyjs_2.1.0 sandwich_3.0-2
[31] labeling_0.4.2 sass_0.4.1 mvtnorm_1.1-3
[34] spatstat.data_3.0-0 brms_2.19.0 proxy_0.4-27
[37] StanHeaders_2.26.13 DHARMa_0.4.6 MuMIn_1.47.5
[40] sessioninfo_1.2.2 readxl_1.4.0 rstudioapi_0.13
[43] circular_0.5-0 visNetwork_2.1.0 generics_0.1.2
[46] gtools_3.9.2.2 crosstalk_1.2.0 vroom_1.5.7
[49] car_3.1-2 distributional_0.3.0 inline_0.3.19
[52] loo_2.5.1 Matrix_1.5-4.1 fansi_1.0.3
[55] abind_1.4-5 lifecycle_1.0.3 multcomp_1.4-25
[58] yaml_2.3.5 carData_3.0-5 grid_4.3.1
[61] qgam_1.3.4 promises_1.2.0.1 crayon_1.5.1
[64] miniUI_0.1.1.1 lattice_0.21-8 haven_2.5.0
[67] cowplot_1.1.1 pillar_1.9.0 knitr_1.43
[70] boot_1.3-28.1 estimability_1.3 shinystan_2.6.0
[73] codetools_0.2-19 imagerExtra_2.0.0 glue_1.6.2
[76] RcppArmadillo_0.11.2.0.0 V8_4.2.0 remotes_2.4.2
[79] vctrs_0.6.3 png_0.1-7 Rdpack_2.3.1
[82] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
[85] cachem_1.0.6 xfun_0.39 rbibutils_2.2.8
[88] mime_0.12 coda_0.19-4 shinythemes_1.2.0
[91] units_0.8-0 DiagrammeR_1.0.9 ellipsis_0.3.2
[94] fitdistrplus_1.1-11 TH.data_1.1-1 gap_1.2.3-6
[97] nlme_3.1-162 CircStats_0.2-6 xts_0.13.1
[100] usethis_2.1.6 bit64_4.0.5 threejs_0.3.3
[103] rstan_2.26.13 rprojroot_2.0.3 tensorA_0.36.2
[106] bslib_0.3.1 TMB_1.9.4 KernSmooth_2.23-20
[109] colorspace_2.0-3 DBI_1.1.3 tidyselect_1.2.0
[112] processx_3.6.1 emmeans_1.7.5 curl_4.3.2
[115] bit_4.0.4 compiler_4.3.1 amt_0.2.1.0
[118] rvest_1.0.2 arrayhelpers_1.1-0 see_0.7.1
[121] xml2_1.3.3 desc_1.4.1 colourpicker_1.1.1
[124] bayestestR_0.13.1 posterior_1.2.2 dygraphs_1.1.1.6
[127] checkmate_2.1.0 scales_1.2.1 classInt_0.4-7
[130] callr_3.7.0 tiff_0.1-11 digest_0.6.29
[133] fftwtools_0.9-11 spatstat.utils_3.0-1 minqa_1.2.4
[136] rmarkdown_2.24 base64enc_0.1-3 htmltools_0.5.2
[139] moveHMM_1.9 pkgconfig_2.0.3 jpeg_0.1-9
[142] lme4_1.1-29 highr_0.9 dbplyr_2.2.0
[145] fastmap_1.1.0 rlang_1.1.1 htmlwidgets_1.5.4
[148] shiny_1.7.1 farver_2.1.0 jquerylib_0.1.4
[151] svUnit_1.0.6 zoo_1.8-12 jsonlite_1.8.8
[154] bayesplot_1.10.0 geosphere_1.5-18 munsell_0.5.0
[157] Rcpp_1.0.8.3 stringi_1.7.12 MASS_7.3-60
[160] plyr_1.8.7 pkgbuild_1.3.1 parallel_4.3.1
[163] sjmisc_2.8.9 deldir_1.0-6 ggeffects_1.2.3
[166] splines_4.3.1 hms_1.1.1 sjstats_0.18.2
[169] ps_1.7.1 igraph_1.3.2 spatstat.geom_3.0-3
[172] markdown_1.1 ggsignif_0.6.3 reshape2_1.4.4
[175] stats4_4.3.1 pkgload_1.3.2 rstantools_2.2.0
[178] reprex_2.0.1 evaluate_0.15 RcppParallel_5.1.5
[181] modelr_0.1.8 bmp_0.3 nloptr_2.0.3
[184] tzdb_0.3.0 httpuv_1.6.5 RgoogleMaps_1.4.5.3
[187] polyclip_1.10-0 broom_1.0.5 gap.datasets_0.0.5
[190] xtable_1.8-4 e1071_1.7-11 rstatix_0.7.2
[193] later_1.3.0 class_7.3-22 memoise_2.0.1
[196] ggmap_3.0.2 bridgesampling_1.1-2